home *** CD-ROM | disk | FTP | other *** search
/ Precision Software Appli…tions Silver Collection 4 / Precision Software Applications Silver Collection Volume 4 (1993).iso / stats / chadyn.exe / YTOOLS.C < prev    next >
Text File  |  1988-12-08  |  20KB  |  765 lines
  1. /******************* (C) 1986,7,8 by JAMES A. YORKE **************************/
  2. /******************************* YTOOLS.C ************************************/
  3.  
  4. /* This file Ytools.C contains routines that are general tools
  5.      that do not require detailed knowledge of 
  6.     the differential equation or map under investigation
  7.             James Yorke               */
  8.  
  9. /*     Routines in YTOOLS.c
  10. These are routines that are essentially independent of the map being studied.
  11.  
  12. beep                
  13. ChooseStorageVec(Character)    Character == 0 means StorVec is not provided
  14. distance(y0,y1,dimen)
  15. error_check             
  16. get_set(output,mapnum)
  17.     int mapnum;
  18.     FILE *output;
  19. iterate_map           
  20. moduloAB(Y,A,B)
  21. MoveVec() depends on level
  22. null() {  }       dummy 
  23. printWhere()
  24. print_text           
  25. print_y()
  26. random
  27. ret_erase_line(i) gives i copies of return(plus line feed) + erase_line() 
  28. rungekutta           
  29. ScanfForCoefs(A,B,pVec)
  30. setequal(J,K,dimen) use store instead of this
  31. start        clears a line & sets dim = 0
  32. store(y0,y1,dimen) sets y0 equal to y1
  33. y_init_init */
  34.  
  35.  
  36.  
  37.  
  38.  
  39.  
  40. /* The core(i.e. pic[]) and screen operations are independent   */
  41. #include "yinclud.h"
  42.  
  43. static long     IterateCounter = 0;
  44. static double   frac =.31415926;/* seed for random() */
  45.  
  46.  
  47. /****************************SUBROUTINES********************************/
  48.  
  49. FILE    *printWhere() 
  50. {
  51.     if(level < PROCESS || prnFlag == OFF) {
  52.         return(StOutPut);
  53.     } else
  54.         return(stdprn);
  55. }
  56. beep() 
  57. {            /* beep() beeps */
  58.     putchar('\007');    /* hopefully goes to screen */
  59. }
  60.  
  61. int     ChooseStorageVec(Character)
  62.                 /* This routine takes a character and puts the
  63.                    address of the corresponding storage vector
  64.                    in ystore_in; it accepts only the characters
  65.                    a, b, ...,e, 0, 1, 2, ... NUM_Y -1, which
  66.                    correspond to storage vectors; if Character
  67.                    == 0 then it gets the character from the
  68.                    keyboard; StorChar is also used by
  69.                    interrupt()  */
  70. char    Character;
  71. {
  72.     int     outcome;
  73.  
  74.     if(Character == 0) {
  75.         /*beep();*/
  76.         if(SCREEN)  {
  77.             if(level >= PROCESS) {
  78.                 PRINT
  79.                     "hit key 0,1,.. < %d OR 'a','b'...'e' indicating which storage vector to use: "
  80.                     ,NUM_Y);
  81. #ifdef X11
  82.                 StorChar = awaitkey(0);
  83. #else
  84.                 while((StorChar = keycheck()) == 0);
  85.                 /* keep trying until a character is entered */
  86. #endif
  87.             }
  88.         }
  89.         if(level == 2) {
  90.             if(SCREEN) {
  91.                 PRINT
  92. "ENTER 0,1,.. < %d OR 'a','b'...'e' THEN <return>                           \n"
  93.                     ,NUM_Y);
  94.                 PRINT
  95. "indicating which storage vector to use:  \n");
  96.             }
  97. #ifdef X11
  98.             if(input == StInput 
  99.                && sscanf(xwfgets(), "%c", &StorChar) != 1
  100.                || input != StInput && 
  101.                    fscanf(input, " %c", &StorChar) != 1)
  102. #else
  103.             if((SCREEN && abortEnter() == YES) || 
  104.                 fscanf(input, " %c", &StorChar) != 1)
  105. #endif
  106.                 return(0);
  107.         }
  108.     } else {
  109.         StorChar = Character;
  110.     }
  111.  
  112.     if(SCREEN)
  113.         PRINT "         \r %c\n", StorChar);
  114.     outcome = 0;
  115.     if((StorChar >= '0' && StorChar < '0' + NUM_Y)
  116.             || (StorChar >= 'a' && StorChar <= 'e'))
  117.         outcome = 1;    /* i.e., character is acceptable */
  118.     else {
  119.         erase_line();
  120.         PRINT " %c is Not acceptable\n", StorChar);
  121.     }
  122.  
  123.  /* NOW DETERMINE WHICH STORAGE VECTOR IS INVOLVED */
  124.     if(StorChar == '0')
  125.         ystore_in = y;
  126.     else if(StorChar >= '1' && StorChar < '0' + NUM_Y)
  127.         ystore_in = y + eqn0 + eqn1 * (StorChar - '1');
  128.                 /* the zeroth vector is allocated eqn0 space,
  129.                    the rest eqns */
  130.     else if(StorChar == 'a')
  131.         ystore_in = ya;
  132.     else if(StorChar == 'b')
  133.         ystore_in = yb;
  134.     else if(StorChar == 'c')
  135.         ystore_in = yc;
  136.     else if(StorChar == 'd')
  137.         ystore_in = yd;
  138.     else if(StorChar == 'e')
  139.         ystore_in = ye;
  140.     return(outcome);
  141. }
  142.  
  143. #ifdef TESTING            /* not defined */
  144. double  dabs(doub)        /* abs() takes int arguments and returns an int
  145.                    */
  146. double  doub;
  147. {
  148.     if(doub < 0)
  149.         doub = -doub;
  150.     return(doub);
  151. }
  152. #endif                /* TESTING */
  153.  
  154. /* yp0 and yp1 are pointers; the distance between
  155.  * the vectors(with dimen coordinates)
  156.  * starting at addresses yp1 and y2 is computed
  157.  * and returned.
  158.  */
  159. double  distance(yp0, yp1, dimen)
  160. double *yp0,
  161.        *yp1;
  162. int     dimen;
  163. {
  164.     double  distsqrd;
  165.     int     i;
  166.  
  167.     distsqrd = 0;
  168.     for(i = 0; i < dimen; i++)
  169.         distsqrd = distsqrd + (yp0[i] - yp1[i]) * (yp0[i] - yp1[i]);
  170.     return(sqrt(distsqrd));
  171. }
  172.  
  173.  
  174. error_check() {        /* This routine is a substitute for rungekutta;
  175.                    it compares results of a dif eq step with
  176.                    two half size steps... The distance between
  177.                    them(=square root of the sum of the squares
  178.                    of the differences of the y[i] values) is
  179.                    the error.  This error is the error for one
  180.                    step only and must be recomputed many times
  181.                    to get an overview of the error size along a
  182.                    trajectory; the final value of y[] is the
  183.                    result of the two half-sized steps; the
  184.                    results are printed on the screen if
  185.                    "printer" > 2 */
  186.     double  y_step_dist,
  187.             distance();
  188.     double  y_one[EQTNS0],
  189.             y_two[EQTNS0];
  190.  
  191.  /* FIRST THE SOLVER IS APPLIED USING THE CURRENT VALUE OF STEP, AFTER SAVING
  192.     THE CURRENT y[] VALUE */
  193.  /* eqns >= dim but this causes no problems since vectors are initialized */
  194.     store(y_one, y, eqn0);    /* store vector y[] in vector y_one[], both
  195.                    having dimension eqns */
  196.     (*IterIteratee) ();    /* this usually equals *map which may equal
  197.                    rungekutta */
  198.  
  199.     y_step_dist = distance(y, y_one, eqn0);
  200.                 /*  this distance between y[] and y_one[] will
  201.                    be used to measure how big the step is in
  202.                    phase space using all variables including
  203.                    any Lyapunov vectors that are being used */
  204.  
  205.     store(y_two, y, eqn0);    /* store vector y[] after one step in vector
  206.                    y_two[], both having dimension eqn0 */
  207.     store(y, y_one, eqn0);    /* store vector y_one[] back in vector y[] */
  208.  
  209.  
  210.  /* NEXT A TWO ITERATES ARE TAKEN WITH HALF AS BIG A VALUE OF STEP */
  211.     step /= 2;
  212.     halfstep /= 2;
  213.     sixth_step /= 2;
  214.     (*IterIteratee) ();
  215.     (*IterIteratee) ();    /* this usually equals *map which may equal
  216.                    rungekutta */
  217.     step *= 2;
  218.     halfstep *= 2;
  219.     sixth_step *= 2;
  220.  
  221.     error = distance(y, y_two, eqn0);
  222.     if(error > max_error)
  223.         max_error = error;
  224.     if(printer > 2 || ODEcheck > 1) {
  225.                 /* ODEcheck is set > 1 by interrupt 'e' but is
  226.                    decreased below in this routine */
  227.         if(num_lyap <= 0)
  228.             scr_rowcol(0, 0);
  229.                 /* move cursor to top left; if num_lyap > 0,
  230.                    the lyap exponents might be printed on top
  231.                    the errors if we returned the cursor */
  232.  
  233.         if(y_step_dist > 0)/* checking so we don't divide by zero */
  234.             PRINT
  235.                 "error= %6.6le;error/step= %6.6le; max error so far = %6.6le    \r"
  236.                 ,error, error / y_step_dist, max_error);
  237.     /* note \r ends line and returns carriage without advancing a line */
  238.     }
  239.     if(ODEcheck > 1)    /* interrupt 'e' increases ODEcheck by 2 so
  240.                    that this routine will be called once; */
  241.         ODEcheck = ODEcheck - 2;
  242. }
  243.  
  244.  
  245. get_set(output, mapname)
  246. char   *mapname;
  247. FILE * output;
  248. {
  249.     if(strcmp(mapname, "all") == 0)
  250.         fputs(setup, output);
  251. }
  252.  
  253. iterate_map() {        /* carries out one iterate(possibly running
  254.                    error check every 5 dots for dif eqn); after
  255.                    iteration it calls(*modPointer)() to see if
  256.                    some value should be taken mod something;
  257.                    this cannot be done by the differential
  258.                    equation solver; for non differential
  259.                    equation maps, the variables are usually
  260.                    calculated modulo whatever directly within
  261.                    the map; finally routine checks to see if
  262.                    y[] should be printed out */
  263.     int     n;
  264.     int     OneFifthIC,
  265.             ICmod5;        /* for checking if IC divisible by 5   */
  266.  
  267.  
  268.     if(num_lyap > 0 && dot >= 0)/* dot = 0 is the case where lyapunov is
  269.                    initialized */
  270.         if(InNewton == 0)/* = 0 if newton is not being called */
  271.             lyapunov();
  272.  
  273.  
  274.     for(n = 0; n < its_per_plot; n++) {
  275.                 /* permits several iterates per point plotted 
  276.                 */
  277.         if(step == -9999.)
  278.             (*IterIteratee) ();
  279.         else {
  280.  
  281.             IterateCounter++;
  282.             OneFifthIC = IterateCounter / 5;
  283.                 /* integer arithmetic */
  284.             ICmod5 = IterateCounter - 5 * OneFifthIC;
  285.  
  286.         /* ODEcheck == 1 is turned on by 'E' in YINTRPT.C */
  287.             if((ODEcheck >= 1) && (ICmod5 == 0))
  288.                 error_check();
  289.                 /* Happens if dot is divisible by 5 */
  290.             else
  291.                 (*IterIteratee) ();
  292.  
  293.             (*modPointer) ();
  294.                 /* this = null() for non differential equations
  295.                    and for many differential equations */
  296.         }
  297.     }
  298. }
  299.  
  300. double  moduloAB(Y, A, B)    /* returns a value that is Y shifted by a
  301.                    multiple of B-A; it is shifted so that the
  302.                    returned value lies between A and B; note A
  303.                    and B and Y must be "double" reals; if the
  304.                    operation changes Y, then modFlag is set =
  305.                    YES = 1; it should be turned off elsewhere;
  306.                    it is not clear to me what happens when a
  307.                    negative double is truncated to an integer
  308.                    so this routine shifts it to be positive
  309.                    before truncating and then shifts back */
  310. double  Y,
  311.         A,
  312.         B;
  313. {
  314.     double  scaledY;
  315.     int     int_part;
  316.     double  shift;
  317.     scaledY = (Y - A) / (B - A);
  318.     if(scaledY >= 0 && scaledY < 1)
  319.         return(Y);
  320.     modFlag = YES;        /* this means that some number is being changed
  321.                    by this routine */
  322.     shift = 10.0;
  323.     while(shift + scaledY < 0)
  324.         shift *= 10.0;
  325.     int_part = shift + scaledY;
  326.                 /* we only take the integer part of a positive
  327.                    number; */
  328.  
  329.     return(A + (B - A) * (scaledY + shift - int_part));
  330. }
  331.  
  332. MoveVec() {            /* asks for source vector and then the target
  333.                    vector, each designated by one character;
  334.                    then the target is set equal to the value in
  335.                    the source. */
  336.     char    Sto;
  337.     int     acceptable;
  338.     double *ToBeMoved;
  339.  
  340.     if(level >= PROCESS) {    /* yintrpt case */
  341.         scr_rowcol(0, 0);/* move cursor to upper left, a routine in
  342.                    desmet's pcio.a */
  343.     }
  344.     else
  345.         PRINT "\n");
  346.  
  347.  
  348.     erase_line();
  349.     PRINT
  350.         "Hit digit or letter of vector whose coordinates will change(Use 0 for y[])\n"
  351.         );
  352.     acceptable = ChooseStorageVec(0);
  353.                 /* Character == 0 means StorVec is not provided
  354.                    */
  355.     Sto = StorChar;
  356.     if(acceptable == 1) {
  357.         if(level >= PROCESS) {/* yintrpt case */
  358.             scr_rowcol(0, 0);
  359.                 /* move cursor to upper left, a routine in
  360.                    desmet's pcio.a */
  361.             erase_line();
  362.             PRINT "\n");
  363.             erase_line();
  364.             PRINT "\n");
  365.             erase_line();
  366.             scr_rowcol(1, 0);
  367.         }
  368.         else
  369.             PRINT "\n");
  370.  
  371.         ToBeMoved = ystore_in;/* pointers */
  372.         if(input == StInput) {
  373.             erase_line();
  374.             PRINT "y%c. Fine. \n", StorChar);
  375.             PRINT
  376.                 "Now hit the digit or letter for vec whose coordinates will be copied to y%c\n"
  377.                 ,StorChar);
  378.         }
  379.         acceptable = ChooseStorageVec(0);
  380.         if(acceptable == YES) {
  381.             store(ToBeMoved, ystore_in, eqn1);
  382.             if(level >= PROCESS) {/* yintrpt case */
  383.                 scr_rowcol(0, 0);
  384.                 /* move cursor to upper left, a routine in
  385.                    desmet's pcio.a */
  386.                 erase_line();
  387.                 PRINT "\n");
  388.                 erase_line();
  389.                 PRINT "\n");
  390.                 erase_line();
  391.                 PRINT "\n");
  392.                 erase_line();
  393.                 scr_rowcol(1, 0);
  394.             }
  395.             PRINT "\n");
  396.             if(input == StInput)
  397.                 PRINT
  398.                     " y%c[] has been copied to y%c[]        \n\n", StorChar, Sto);
  399.             ystore_in = ToBeMoved;
  400.             return;
  401.         }
  402.     }
  403.     PRINT "Not acceptable; start over to try again \n");
  404. }
  405.  
  406.  
  407. null() {
  408. }                /* dummy */
  409.  
  410.  
  411. print_text(output)        /* for sending output of data to file output;
  412.                    in particular it can go to printer LPT1
  413.                    which is done by calling
  414.                    print_text(StPrint); at least this works
  415.                    with IBMPCs; the amount of data printed
  416.                    depends on "printer" or to the screen */
  417. FILE * output;
  418. {
  419.     double  x,
  420.             X,
  421.             yval,
  422.             Y;
  423.  
  424.     if(printer > 0)
  425.         if(fprintf(output, " ") == -1) {/* check for error */
  426.             printf(
  427.                     "The program cannot access the output file. Perhaps the printer is off \n");
  428.             return;
  429.         }
  430.  
  431.     if(bifTextFlag == OFF)    /* i.e. not using command BIF */
  432.         if(printer == 1) {
  433.             if(rho != -9999.)
  434.                 fprintf(output, "rho = %15.10lf        ", rho);
  435.             if(output == StOutPut)
  436.                 fprintf(output, "Last  dot was %ld. ", dot);
  437.             fprintf(output, "\n");
  438.         }
  439.     if(printer > 1) {
  440.         if(bifTextFlag == OFF)/* i.e. not using command BIF */
  441.             if(picNameFlag == YES)
  442.                 /* means DiskFileName has been accessed, either
  443.                    storing or reading */
  444.                 fprintf(output,
  445.                         "Associated Disk file name: \"%s\" \n", DiskFileName);
  446.         map_menu(output, MapName);
  447.         if(bifTextFlag == OFF)/* i.e. not using command BIF */
  448.             if(rho != -9999.)
  449.                 fprintf(output, "rho = %lf     ", rho);
  450.         if(step != -9999)
  451.             fprintf(output, "step = %lf\n", step);
  452.         if(C1 != -9999)
  453.             fprintf(output, "C1 = %lf  ", C1);
  454.         if(C2 != -9999)
  455.             fprintf(output, "C2 = %lf  ", C2);
  456.         if(C3 != -9999)
  457.             fprintf(output, "C3 = %lf  ", C3);
  458.         if(C4 != -9999)
  459.             fprintf(output, "C4 = %lf  ", C4);
  460.         if(C5 != -9999)
  461.             fprintf(output, "C5 = %lf  ", C5);
  462.         if(C6 != -9999)
  463.             fprintf(output, "C6 = %lf  ", C6);
  464.         if(C7 != -9999)
  465.             fprintf(output, "C7 = %lf  ", C7);
  466.         if(C8 != -9999.)
  467.             fprintf(output, "C8 = %lf  ", C8);
  468.         if(beta != -9999.)
  469.             fprintf(output, "  beta  = %lf ;", beta);
  470.         if(sigma != -9999.)
  471.             fprintf(output, "  sigma = %lf ;", sigma);
  472.         fprintf(output, "\n");
  473.  
  474.         x = X_Lo[ScrnSec];
  475.         X = X_Up[ScrnSec];
  476.         yval = Y_Lo[ScrnSec];
  477.         Y = Y_Up[ScrnSec];
  478.  
  479.         if(ScrnSec != 0)
  480.             fprintf(output, "Section F%d X coord: from %lf to %lf; "
  481.                     ,ScrnSec, x, X);
  482.         else
  483.             fprintf(output, "X coord: from %lf to %lf; ", x, X);
  484.         if(bifTextFlag == OFF)/* i.e. not using command BIF */
  485.             fprintf(output, "Y coord: %lf to %lf\n", yval, Y);
  486.         if(preiter != 0)
  487.             fprintf(output, "preiterates before plotting = %ld\n"
  488.                     ,preiter);
  489.         if(bifTextFlag == OFF)/* i.e. not using command BIF */
  490.             fprintf(output, "Last  dot was %ld.   ", dot);
  491.  
  492.         if(bifTextFlag == OFF)/* i.e. not using command BIF */
  493.             if(dim > 2 || X_coord != 0 || Y_coord != 1)
  494.                 fprintf(output,
  495.                         "Coordinate numbers(0 to %d) are: %d and %d"
  496.                         ,dim - 1, X_coord, Y_coord);
  497.         fprintf(output, "\n");
  498.         if(lyaptime > 0.&& num_lyap > 0)
  499.             print_lyap(output);
  500.         if(output != StOutPut)
  501.             fprintf(output, "\n\n");
  502.     }            /* end > 1 case */
  503. }
  504.  
  505.  
  506. print_y() {            /* used with newton method, called from
  507.                    YCASES2.C */
  508.     int     ii;
  509.  
  510.     for(ii = 0; ii < lyapzero; ii++) {
  511.         erase_line();
  512.         PRINT "y[%d] = %15.12lf      \n", ii, y[ii]);
  513.     }
  514.     if(vec_dim == 2 && lyapzero > 2) {
  515.         erase_line();
  516.         PRINT
  517.             "            The phase space variables are y[%d] and y[%d]\n"
  518.             ,zeroth, zeroth + 1);
  519.     }
  520. }
  521.  
  522. int     random997 () {        /* Returns an integer between 0 and 997 using
  523.                    external frac */
  524.     double  rdm;
  525.     int     R;
  526.  
  527.     rdm = 997 * frac;
  528.     R = rdm;
  529.     frac = rdm - R;
  530.     return(R);
  531. }
  532.  
  533. euler() {            /* this Euler routine is for a single step of
  534.                    the differential equation; DE is a pointer
  535.                    to a function  */
  536.     int     i;
  537.     double  k[EQTNS0];
  538.  
  539.     (*DE) (k, y);
  540.     for(i = 0; i < dim; i++)
  541.         y[i] += step * k[i];
  542.     t = t + step;
  543. }
  544.  
  545. ret_erase_line(i)        /* gives i copies of return(plus line feed) +
  546.                    erase_line() */
  547. int     i;
  548. {
  549.     int     n;
  550.  
  551.     for(n = 0; n < i; n++) {
  552.         PRINT "\n");
  553.         erase_line();
  554.     }
  555. }
  556.  
  557. rungekutta() {            /* this fourth order Runge_Kutta routine is for
  558.                    a single step of the differential equation;
  559.                    DE is a pointer to a function  */
  560.     int     i;
  561.     double  k1[EQTNS0],
  562.             k2[EQTNS0],
  563.             k3[EQTNS0],
  564.             k4[EQTNS0];
  565.     double  temp[EQTNS0];
  566.     (*DE) (k1, y);
  567.     for(i = 0; i < dim; i++)
  568.         temp[i] = y[i] + halfstep * k1[i];
  569.     (*DE) (k2, temp);
  570.  
  571.     for(i = 0; i < dim; i++)
  572.         temp[i] = y[i] + halfstep * k2[i];
  573.     (*DE) (k3, temp);
  574.  
  575.     for(i = 0; i < dim; i++)
  576.         temp[i] = y[i] + step * k3[i];
  577.     (*DE) (k4, temp);
  578.     for(i = 0; i < dim; i++)
  579.         y[i] += sixth_step * (k1[i] + 2.* k2[i] + 2.* k3[i] + k4[i]);
  580.     t = t + step;
  581. }
  582.  
  583.  
  584. int     ScanfForCoefs(A, B, pVec)
  585.                 /* the coordinates acceptable are >= A and < B;
  586.                    the routine asks for a coordinate first; if
  587.                    it is acceptable, it asks for a double
  588.                    precision value to plug into the
  589.                    corresponding coordinate; if successful, and
  590.                    coord < B-1, a 1 is returned so inputting
  591.                    can continue; otherwise 0 is returned; if
  592.                    the routine is not successful in getting
  593.                    both a coordinate and a double, pVec is not
  594.                    changed */
  595. int     A,
  596.         B;
  597. double *pVec;
  598. {
  599.     int     coord,
  600.             i;
  601.     double  DOUB;
  602.     int     ret;
  603.  
  604.     if(input == StInput)
  605.         for(i = A; i < B; i++)
  606.             PRINT "coordinate %d = %10.10lf \n", i, pVec[i]);
  607.  
  608.     if(input == StInput)
  609.         puts("To change a coordinate, input coordinate number and RETURN\n");
  610.     if(input == StInput)
  611.         puts("To LEAVE values UNCHANGED, just hit RETURN\n");
  612.     coord = -1;
  613.     DOUB = -19999.;
  614.     ret = 0;
  615. #ifdef X11
  616.     if(input == StInput && sscanf(xwfgets(), "%d", &coord) != 1
  617.       || input != StInput && fscanf(input, "%d", &coord) != 1)
  618.          return(0);
  619. #else
  620.     if (SCREEN && abortEnter() == YES || fscanf(input," %d",&coord) != 1)
  621.         return(0);
  622. #endif
  623.     if(coord < B && coord >= A) {
  624.         if(input == StInput)
  625.             PRINT
  626.                 "Enter new value of coordinate # %d, then RETURN\n", coord);
  627.  
  628. #ifdef X11
  629.         if(input == StInput && sscanf(xwfgets(), "%lf", &DOUB) != 1
  630.           || input != StInput && fscanf(input, "%lf", &DOUB) != 1)
  631.             return(0);
  632. #else
  633.         if(SCREEN && abortEnter() == YES ||
  634.             fscanf(input, "%lf", &DOUB) != 1)
  635.             return(0);
  636. #endif
  637.         if(DOUB != -19999) {
  638.             pVec[coord] = DOUB;
  639.             ret = 1;/* means coordinate was satisfactorily entered 
  640.                 */
  641.             if(coord == B - 1)
  642.                 ret = 0;/* B-1 = last coord */
  643.         }
  644.     }
  645.     else if(coord != -1)
  646.         puts("Illegal coordinate. Not accepted. \n");
  647.     if(ret == 0)
  648.         if(input == StInput)
  649.             for(i = A; i < B; i++)
  650.                 PRINT "coordinate %d = %lf\n", i, pVec[i]);
  651.     return(ret);
  652. }
  653.  
  654.  
  655. setequal(J, K, dimen)
  656. int     J,
  657.         K,
  658.         dimen;
  659. {
  660.     int     i,
  661.             j,
  662.             k;
  663.  
  664.     if(J == 0)
  665.         j = 0;
  666.     else
  667.         j = eqn0 + eqn1 * (J - 1);
  668.  
  669.     if(K == 0)
  670.         k = 0;
  671.     else
  672.         k = eqn0 + eqn1 * (K - 1);
  673.     for(i = 0; i < dimen; i++)
  674.         y[j + i] = y[k + i];
  675. }
  676.  
  677.  
  678.  
  679. store(yp0, yp1, dimen)        /* yp0 and yp1 are pointers; the vector starting
  680.                    at yp1 is stored in the vector beginning at
  681.                    yp0; dimen is the number of coordinates to be
  682.                    transfered */
  683. double *yp0,
  684.        *yp1;
  685. int     dimen;
  686. {
  687.     int     i;
  688.     for(i = 0; i < dimen; i++)
  689.         yp0[i] = yp1[i];
  690. }
  691.  
  692. y_init_init() {        /* this routine sets the initalization vector;
  693.                    it is used before the default parameters of
  694.                    a map or differential equation are set;
  695.                    hence other values can be set when the
  696.                    defaults are set */
  697.     int     i,
  698.             k;
  699.     for(i = 0; i < eqn1; i++) {
  700.         y[eqn0 + i] = 0.;
  701.         for(k = 2; k < NUM_Y; k++)
  702.             y[eqn0 + (k - 1) * eqn1 + i] = -9999.;
  703.     }
  704.     y[eqn0 + 1] = 1.;
  705. }
  706.  
  707.  
  708.  
  709. double  Fn3 () {        /* this is designed for the Lorenz system; the
  710.                    pair of stationary solutions that the
  711.                    solution oscillates around has z(= y[2]) =
  712.                    rho - 1 */
  713.             return(y[2] - rho + 1);
  714. }
  715.  
  716. double  FnLinY() {
  717.     int     J;
  718.     double  F;
  719.     F = 0;
  720.     for(J = 0; J < dim; J++)
  721.         F = F + PoincareCoef[J] * y[J];
  722.     return(F);
  723. }
  724. /* insert and test later
  725. double FnDeriv()... this is a linear function of the derivative with the 
  726.         same coefficients as FnLinY() 
  727. {
  728.     int J;
  729.     double F;
  730.     double k1[EQTNS0];
  731.     (*DE)(k1,y);
  732.     F = 0;
  733.     for(J = 0; J < dim; J++)
  734.         F = F + PoincareCoef[J]*y[J];
  735.     return(F);
  736. }
  737. */
  738.  
  739. SetPoincareCoefs() {
  740.     int     J;
  741.     if(input == StInput)
  742.         PRINT
  743.             "set Poincare cross section coeficients for J >= 0, J < %d; currently \n", dim);
  744.     if(input == StInput)
  745.         for(J = 0; J < dim; J++)
  746.             PRINT
  747.                 " PoincareCoef[%d] = %lf;         \n", J, PoincareCoef[J]);
  748.     while(ScanfForCoefs(0, dim, PoincareCoef)) {;
  749.     }
  750.     return;
  751. }
  752. #ifndef UNIX
  753. pause(sec)            /* pause for this number of seconds */
  754. double  sec;            /* we use sleep() on Unix systems */
  755. {
  756.     int     timevector[4];
  757.     double  timeofday(), time0;
  758.  
  759.     time0 = timeofday(timevector);
  760.  
  761.     while(timeofday(timevector) - time0 < sec);
  762.     return;
  763. }
  764. #endif    /* UNIX */
  765.